#coding=utf-8
'''
Created on 2011-3-18

@author: zarra
'''
import miniDicom as dicom
import struct
import numpy as np
from pylab import *
import matplotlib.cm as cm
import matplotlib as mpl
import scipy.ndimage.morphology as m
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False


dicom_image = dicom.read_file('../data/3tm-1.dcm')

height = dicom_image.Rows
width  = dicom_image.Columns
size = height*width
bits_allocated = dicom_image.BitsAllocated
bits_stored    = dicom_image.BitsStored

#获得dicom图像矩阵
data = dicom_image.PixelData
data = np.array(struct.unpack('<%dH'%size,data),dtype=np.int)
data = data.reshape(height,width)


  
#以20为阈值进行二值化
binarization = np.frompyfunc( lambda  value:  1 if value > 20 else 0,1,1)
b_data = binarization(data)
b_erode= np.copy(b_data).astype(np.int)

E1=9 
E2=9

#用于腐蚀的结构体,9x9的方阵
struct1=np.ones((E1,E1)).astype(np.int)
#用于膨胀的结构体,9x9的方阵
struct2=np.ones((E2,E2)).astype(np.int)

#三次腐蚀
b_erode=m.binary_erosion(b_erode, struct1).astype(np.int)
b_erode=m.binary_erosion(b_erode, struct1).astype(np.int)
b_erode=m.binary_erosion(b_erode, struct1).astype(np.int)


#三次膨胀
b_erode=m.binary_dilation(b_erode, struct2).astype(np.int)
b_erode=m.binary_dilation(b_erode, struct2).astype(np.int)
b_erode=m.binary_dilation(b_erode, struct2).astype(np.int)
     


    
    
smooth = np.copy(data)

w=0
h=0
MASK=3
#根据二值数据重建图像
while h <height:
    while w<width: 
        smooth[h,w] = smooth[h,w] if b_erode[h,w] else 0
        #print v
        w+=1
    h+=1
    w=0    


edge=np.zeros((height,width,3))
w=0
h=0
MAX = np.max(data)
#Roberts边缘检测算法 |z_9-z_5|+|z_8-z_6|
while h <height-MASK:
    while w<width-MASK: 
        M = b_erode[h:h+MASK,w:w+MASK]
        v =np.abs(M[2,2]-M[1,1])+np.abs(M[2,1]-M[1,2])
        vi = data[h,w]/float(MAX)
        edge[h,w,:]=vi
        if v >0:
            edge[h,w,0] =1.0
            edge[h,w,1] =0.0
            edge[h,w,2] =0.0
        w+=1
    h+=1
    w=0  

figure(figsize=(8,4))
subplot(221)
imshow(data,cmap=cm.Greys_r)
title(u"原图 ")


subplot(222)
imshow(smooth,cmap=cm.Greys_r)
title(u"去边缘")

subplot(223)
imshow(data-smooth,cmap=cm.Greys_r)
title(u"差图像")

subplot(224)
imshow(edge)
title(u"边缘")
show()